function lambda = comp_eigs_adj(k, p, h, X0)

if nargin == 3
    X0 = [0; 0];
end

opts = optimset('tolfun',1e-10, 'tolX',1e-10);
X = fsolve(@(x) fun(x, h, k, p), X0, opts);

fprintf('lambda = %.4e + 1i * %.4e\n', X(1), X(2));

lambda = X(1) + 1i * X(2);

function F = fun(X, h, k, p)

lambda = X(1) + 1i * X(2);

mu = sqrt(-k^2 - lambda);

A = w_coeffs(1, h, k, p);
w = @(z)  (A(1) * z + A(2)) .* cosh(k*z) + (A(3) * z + A(4)) .* sinh(k*z);

v_0 = (1 - p.beta) * (1 - 1 / h);
c = p.beta + v_0 - p.delta * (1-p.beta)*(p.beta + v_0) / 2;
gamma = p.delta * (1 - 2 * c);


% the problem in general
tmp = mu * sin(mu * h) - gamma * cos(mu * h) + ...
    p.delta * (1-p.beta) * (p.beta + v_0) / h^2 * quad(@(z) w(z) .* z .* cos(mu*z), 0, h, 1e-10);

% in the large Ma limit:
% tmp =  quad(@(z) w(z) .* z .* cos(mu*z), 0, h);


F(1) = real(tmp);
F(2) = imag(tmp);